Three-sublattice order in the SU(3) Heisenberg model on the square and triangular lattice 
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We present a numerical study of the SU(3) Heisenberg model of three-flavor fermions on the triangular and 
square lattice by means of the density-matrix renormalization group (DMRG) and infinite projected entangled- 
pair states (iPEPS). For the triangular lattice we confirm that the ground state has a three-sublattice order with 
a finite ordered moment which is compatible with the result from linear flavor wave theory (LFWT). The same 
type of order has recently been predicted also for the square lattice [PRL 105, 265301 (2010)] from LFWT and 
exact diagonalization. However, for this case the ordered moment cannot be computed based on LFWT due to 
divergent fluctuations. Our numerical study clearly supports this three-sublattice order, with an ordered moment 
of m = 0.2 — 0.4 in the thermodynamic limit. 

PACS numbers: 67.85.-d, 71.10.Fd, 75.10.Jm, 02.70.-c 



I. INTRODUCTION 

Recent advances in experiments on cold atomic gases have 
raised interest in systems consisting of several flavors of inter- 
acting fermions, which can be realized for example as differ- 
ent hyperfine states of alkali atoms 1 ,2 or nuclear spin states of 
ytterbium or alkaline-earth atoms. 1,4 A model Hamiltonian 
to describe such systems is the TV-flavor fermionic Hubbard 
model given by 



H.c. 



riiaUip, (1) 



where a, (3 run over the different flavors, (i, j) runs over pairs 
of nearest neighbors on the lattice and i runs over all sites of 
the lattice. 

The two-flavor case corresponds to the spin-^ Hubbard 
model. It is generally accepted that for sufficiently large U 
and at half filling, i.e. when each lattice site is occupied by 
exactly one fermion, the ground state is an antiferromagnetic 
Mott insulator. In experiments on cold atoms, the transition 
to a Mott insulator has recently been observed; v the obser- 
vation of the antiferromagnetic spin order is still open. The 
spin-| Heisenberg model is believed to be a good low-energy 
model for the spin degrees of freedom. 

For the more general case N > 2, it is expected that, at cer- 
tain fillings, Mott insulating states will also emerge. ' How- 
ever, the spin order (or flavor order) in this case is not under- 
stood yet. This has raised interest in generalizations of the 
spin-| Heisenberg model, namely SU(A) Heisenberg mod- 
els. Analogous to the N — 2 case, these are obtained as 
second-order expansion of the above A-flavor fermionic Hub- 
bard model in t/U at a filling such that each site is occupied 
by exactly one particle. The Hamiltonian is 



(2) 



where the first sum runs over pairs of nearest neighbors and 
the second sum over flavors. In this paper, we will focus 
on the case of the SU(3) Heisenberg model, where a, f3 g 
{A,B,C}. 

Note that this model is different from the SU(AT) Heisen- 
berg models studied in Refs. 9-17 where other irreducible 
representations of SU(AT) have been considered, which can 
be labelled by different Young tableaus. The corresponding 
Young tableau of our model has one single box, i.e. the fun- 
damental representation at each site. In Refs. and 10 the 
large-N limit with N/2 particles per site has been studied (a 
Young tableau with N/2 boxes in one column). In Refs. 16 
and 17 the large-N limit with representations with m rows and 
n c columns with fixed N/m and n c are considered (n c = 1 in 
Ref. 16). Another possibility is to use conjugate representa- 
tions on two sublattices, 10-1 ^ which is accessible by Quantum 
Monte Carlo simulations without a sign problem, ~ in con- 
trast to the SU(3) model considered in this paper. 

This SU(3) model is equivalent to the spin-1 bilinear- 
biquadratic model, 



H = ^2 [cos6(Si ■ Sj) + sin 0$ ■ Sj 

(id) 



(3) 



(i,j) a, /3 



with = 7r/4, thus when bilinear and biquadratic terms are 
equal and positive. 

In a pioneering work, Papanicolaou studied the phase di- 
agram of this model on the square lattice as a function of 9 by 
a semiclassical analysis. Using a site-factorized variational 
ansatz (product state) he proposed that the case 8 — tt/4 cor- 
responds to a phase transition from the antiferromagnetically 
ordered phase (adiabatically connected to the purely bilinear 
case) to a "semi-ordered phase" with infinitely many degen- 
erate ground states, including states with 2- or 3-sublattice or- 
der. For recent progress on the nature of the phases in the 
"semi-ordered" region and its vicinity, see Ref. 19. 

The situation changes for the closely related case of the tri- 
angular lattice (equivalent to introducing a coupling on one of 
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FIG. 1. Proposed three-sublattice order for the SU(3) Heisenberg 
model on the square and triangular lattice. The triangular lattice 
is obtained from the square lattice by adding couplings along the 
dashed bonds shown above. Blue boxes indicate the sites that are 
pinned to a specific flavor in our DMRG simulations in order to ex- 
plicitly break SU(3) symmetry. 



the diagonals of each plaquette of the square lattice), where 
a site-factorized ansatz already predicts a three-sublattice or- 
der as depicted in Fig. I. 2 "' 21 This state is stable upon adding 
quantum fluctuations at the level of linear flavor wave theory 
(LFWT), and is supported by exact diagonalization results. 
In contrast, on a zig-zag chain (the one-dimensional analog 
of the triangular lattice) the system undergoes spontaneous 
trimerization. 22 

A recent study based on LFWT indicates that on the square 
lattice a similar type of three-sublattice order is selected by 
quantum fluctuations. This type of order is further supported 
by exact diagonalization revealing a tower of states compati- 
ble with the continuous symmetry breaking of SU(3). 23 The 
ordered moment of the symmetry broken state, however, can- 
not be computed at the level of LFWT because fluctuations are 
divergent. Note that at this point this is an artifact of the linear 
spin wave theory, and it is open whether higher order flavor 
wave corrections would lead to a finite or absent ordered mo- 
ment. The only estimate of the ordered moment so far was 
obtained with exact diagonalization from the real-space cor- 
relation functions of an 18-site cluster, suggesting an ordered 
moment of 60% — 70% of the saturation value, which is ex- 
pected to decrease with system size. To further establish the 
three-sublattice order on the square lattice it is important to 
have an estimate of the ordered moment in the thermodynamic 
limit. 

While on both lattices three-sublattice order has been sug- 
gested, the mechanism how this order is selected is quite dif- 
ferent: on the triangular lattice it is already favored at the 
classical level, and quantum fluctuations only renormalize the 
ordered moment, which is a situation similar to that of the 
SU(2) Heisenberg model on bipartite lattices. In the case of 
the square lattice, on the other hand, the three-sublattice order 
is one among the many degenerate states in the classical limit, 
which is selected by quantum fluctuations. Note that thermal 
fluctuations may select a different order. 

In the previous semiclassical studies quantum fluctuations 
at the level of LFWT have been taken into account, and higher 
order terms have been neglected. Thus, it is still an open 



question if the three-sublattice order is stable upon including 
higher-order quantum fluctuations, or if in this case another 
state is selected. An example of such a scenario has recently 
been observed in the SU(4) Heisenberg model on the square 
lattice, where low-order quantum fluctuations select a pla- 
quette state, but additional higher-order quantum fluctuations 
finally favor a dimerized state. For the SU(3) model, exact re- 
sults on small systems suggest that the order is stable, 21 ' 23 but 
an accurate numerical study for larger systems in two dimen- 
sions is so far missing. 

In this paper we study the stability of the three-sublattice or- 
der of the model (2) on the triangular and square lattice with 
state-of-the-art numerical simulations. We present results for 
finite 2D systems with open boundaries up to a size 8x8 using 
the density matrix renormalization group (DMRG) method, 
and infinite 2D systems with infinite projected entangled-pair 
states (iPEPS). Both methods belong to the class of tensor net- 
work algorithms, enabling to compute ground state properties 
with an accuracy which can be systematically controlled by a 
refinement parameter, called the bond dimension. Both meth- 
ods confirm that the ground state has three-sublattice order 
for both type of lattices, and we provide an estimate of the or- 
dered moment in the thermodynamic limit. Finally we discuss 
an alternative approach based on Schwinger bosons. Unfortu- 
nately, this approach turns out to be unable to describe spon- 
taneous SU(3) symmetry breaking, and, as a consequence, its 
results disagree with those of all other approaches regarding 
the type of ordering and the value of the ordered moment. 

The outline of this paper is as follows. In Sec. II we give 
a short summary of linear flavor wave theory and provide de- 
tails on the DMRG and iPEPS simulations. In Sec. Ill we 
first present the results for the triangular lattice, where the 
three-sublattice order is expected to be more robust than on 
the square lattice, since this order is already favored at the 
classical level. We compare and discuss results for the ener- 
gies and the ordered moment obtained with DMRG, iPEPS 
and LFWT. In Sec. IIIC we provide a similar study for the 
square lattice case, where we find an ordered moment which 
is also finite, but stronger suppressed by quantum fluctuations 
than on the triangular lattice. Finally, Sec. V summarizes our 
results. In Appendix A we report on our attempt to extend the 
Schwinger boson mean-field theory to SU(3). 



II. METHODS 
A. Linear flavor wave theory 

The linear flavor wave theory is the extension of the usual 
SU(2) spin wave theory to SU(N) models. It has been for- 
mulated in Refs. 25 and 18 for the SU(3) case and in Ref. 26 
for the SU(4) case. For completeness, here we give some de- 
tails for the three-sublattice order on the triangular and square 
lattice in the SU(3) Heisenberg model — the cases under 
scrutiny in this paper. We note that for the triangular lattice, 
an analogous calculation has been presented by Tsunetsugu 
and Arikawa. 20 

We begin by extending the Hamitonian (2) to the case 
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where on each site the states belong to the symmetrical ir- 
reducible representation of the SU(3) algebra that can be rep- 
resented by Young-tableaux drawn with M boxes arranged 
horizontally. The SU(3) spin operators in such a symmetrical 
irreducible representation can be expressed as 



S£(/) = b} (l)b a (l), 



(4) 



using Schwinger bosons with 3 flavors, where / is the site in- 
dex and the number of bosons on each site is 



£ 

a£{A,B,C} 



bUi)b a (i) = m, 



(5) 



equal to the number of boxes in the Young tableau. The Sp(l) 
operators satisfy the SU(3) Lie algebra, 



dp i dp> — dp Op, 



na CQ 
dp' O fj 



(6) 



where 6%* is the Kronecker 6 function. For M = 1 the Sg(l) 
operators act on the 3-dimensional, fundamental representa- 
tion | a) (where a — A,B, or C) of the SU(3) algebra as 
S<* \a) = |/3) and Sp*\a') = if a' ^ a, with i being the site 
index. 

The Hamiltonian now can be written as 



H = jJ2Sp(i)SP(3) 

{id) 



(7) 



where the sum is over the nearest neighbor lattice sites, and 
over the repeated a and j3 flavor indices. To draw a parallel to 
the SU(2) case, the Hamiltonian (2) in the fundamental irre- 
ducible representation corresponds to the spin- 1/2 Heisenberg 
model, while the Hamiltonian (7) to the Heisenberg model of 
spins with length S (actually for the SU(2) case S = M/2). 

In the following, we consider an ordered state where the 
spins on the sites I, that belong to sublattice A a , point in the 
direction a. Following the analogy with the spin wave theory 
that is a 1/S expansion, we take the M — > oo limit and do a 
1 /M expansion. Starting from the ordered state we can use 
the following expansion for the Sg(l) operators in the large- 
M limit: 



S%(l)=M-ti a (l), _ 
S%(1) = bf(l)y/M-n a (l)*s VMbf(l), 
S?(l) = y/M-n a (l)b%(l) « yffibftl), 
S%'(l) = bf(l)b a p,(l) 7 

where we have introduced the shorthand notation 



P^a 



(8) 
(9) 
(10) 
(11) 

(12) 



The bp (I) operators with j3 ^ a now correspond to the 
Holstein-Primakoff bosons on sublattice A a , and the super- 
script a keeps track of the sublattice. We replace the expres- 
sions above into Hamiltoniam (7). Expanding in 1/M and 



keeping the quadratic terms only, for the exchange term be- 
tween sites / G A a and I' £ A a , we get 

J2s}(i)s?(i') = M[b a j(i)b«(i)+b a a 'Hnb«'(i>) 

0,7 

+b$(i)b^(i') + b«,(iK'(i%W 

in leading order in M — note that the bosons with flavor dif- 
ferent from the ordered a and a' flavor are missing from the 
bond expression. Assuming a three-sublattice ordered state, 
we define the following Fourier transformation: 

l£A a 

where the summation is over the Aa/3 sites of the A a sub- 
lattice (N\ is the number of lattice sites). The Hamiltonian 
between the sublattices A a and Ap in k-space reads 



Hap — 



zJM 



!,"t 



(15) 



where z is the coordination number of the lattice (z = 4 for 
the square and z = 6 for the triangular lattice). The factor 7k 
reads 



7k = i( e ^+2e-^/ 2 cos^| fc ' 



u 



for the triangular lattice and 
1 , 

7k = x ( 



) 



(16) 



(17) 



for the square lattice, with 7^ = 7 _ k . 

The full Hamiltonian is H = ^2 a< p H a p- It can be diago- 
nalized via a Bogoljubov transformation: 



m 

°a,k 
°/3,-k 



cosh6»(k) sinh(9(k) 
sinh0(k) cosh (9 (k) 



(18) 



with tanh 26(k) = 7 k , leading to 



H = --JMN A + M £ ££^(k) 

kGRBZ a p^ a 



°/3,k°/3,k + 2 



The dispersion of the flavor waves is given by 



(19) 



(20) 



There are 6 degenerate branches in the reduced Brillouin zone, 
which is equivalent to 2 branches in the normal Brillouin 
zone. The dispersion agrees with the result of Tsunetsugu and 
Arikawa for the triangular lattice. For the square lattice, it is 
given in Ref.23. 
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The energy per site due to quantum fluctuations is given by 
the expression 

'w(k)\ 



BZ 



J\M, 



(21) 



where we take into account that there are two modes per lattice 
denotes the average over the Brillouin zone. 



1 BZ 



site. The 

Quantum fluctuations lower the energy from to —0.630 J per 
site for the triangular and to —0.727 J for the square lattice. 
Note that the energy per site of the triangular lattice is higher 
than the one of the square lattice despite the larger coordina- 
tion number of the former lattice. 

The reduction of the ordered moment is calculated as 



(52(0) = M - {^(l)) = M 



1 



1 



' BZ 

(22) 

In the triangular lattice (S%(1)) = M - 0.516, so that the 
on-site moment is reduced from 1 to 0.484. In the square 
lattice, the reduced moment diverges due to the zero line in 
the spectrum. Thus, LFWT is unable to make a prediction for 
the ordered moment. We have tried to use a Schwinger boson 
mean-field theory (SBMFT) to restore a gap along this line 
and remove the divergence (see Appendix A). Unfortunately, 
the SBMFT turned out to be unsatisfactory in several respects, 
and its results regarding the ordered moment are not reliable. 



B. DMRG for finite two-dimensional systems 

1. Setup 

For our DMRG simulations, we map the two-dimensional 
system to a chain following a "TV screen" method (sweeping 
along the vertical direction). We will generally refer to the 
extent in the horizontal direction as length, and in the vertical 
direction as width of the system. We use a single-site opti- 
mization scheme augmented by the improvement suggested in 
Ref. 27. We perform the simulation starting from different ini- 
tial states and increase the bond dimension very quickly with 
the number of sweeps to avoid getting trapped in local min- 
ima. This is particularly important for the case of the square 
lattice, where an insufficient bond dimension may lead to un- 
physical states. Together with the large number of operators, 
this limits the bond dimension that we can reach with our com- 
putational resources to about D ~ 5000 states. Due to the 
huge growth of entanglement with the width of the system, 
this allows us to obtain sufficient accuracy for systems up to 
width 8. 



2. Calculation of the order parameter 

We expect that, in the thermodynamic limit, the SU(3) sym- 
metry is spontaneously broken. If an appropriate basis is cho- 
sen (i.e. after an appropriate SU(3) rotation), one flavor be- 
comes stronger on each site, i.e. 



In this case, we can define the local moment 



( m ) = 7 ) [ max (n a ) - \ 

2 \a=A,B,C 6 



(24) 



which should acquire a finite value in the range (m) £ [0,1]. 

On finite systems, the symmetry is never broken sponta- 
neously and one would conventionally use the relation 



(n a ) 2 = lim ((n a ,in a ,i + 3 C / 

d— r-OO 



) - (n a ,i)(n ati+ 3d)) (25) 



> Tip 



a,/3,j€{A,B,C}. 



(23) 



to extract information about the moments. This however re- 
quires large systems and very accurate estimates for the cor- 
relation functions, which are hard to obtain from a DMRG 
simulation in two dimensions. We therefore follow the pre- 
scription of Refs. and and break SU(3) symmetry ex- 
plicitly by introducing fields at the boundaries of the system. 
The local moments can then be measured locally, preferably 
on sites far away from the pinning fields. The pinning fields 
also fix the direction of the symmetry breaking to be along the 
basis vectors. 

We introduce a column of pinned sites at each end of the 
system, as shown in Fig. 1. We choose the system sizes such 
that the unpinned sites form a square, i.e. the system size in- 
cluding pinned sites is (L + 2) x L, Pinning is done with 
a flavor-specific chemical potential of magnitude 1. In addi- 
tion, such pinning fields reduce the entanglement in the sys- 
tem. Simulations were performed for both open and cylindri- 
cal boundary conditions. 



3. Boundary conditions 

An important question when performing finite-size DMRG 
simulations is the appropriate choice of boundary conditions. 
From an entanglement point of view, open boundary condi- 
tions appear favorable; also, these will have fewer long-range 
operators in the mapping to a chain. From a physical point 
of view, on the other hand, periodic boundary conditions are 
often preferred as they eliminate boundary effects. A compro- 
mise suggested e.g. in Ref. 28 is to use cylindrical boundary 
conditions, which are favorable from an entanglement point 
of view. 

Physically, such boundary conditions are compatible with 
the approach of pinning two columns, which preserves trans- 
lational invariance in the vertical direction. In order not to 
frustrate the three-sublattice order, such boundary conditions 
should only be chosen for systems whose width is a multiple 
of three. For other system sizes, shifted cylindrical boundary 
conditions can be used. For example, for a system of width 5, 
the bottom site of column i must be connected to the top site 
of column i + 1 to obtain a system without additional frustra- 
tion. 

We find numerically that cylindrical boundary conditions 
favor a state that is a product of periodic length-6 chains 
wrapped around the cylinder. A calculation for the same 
model on a periodic chain of length 6 shows that the energy 
per site is very low in this case, making it favorable for small 
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clusters. Such a state shows significantly reduced local mo- 
ments. By choosing open boundary conditions in all direc- 
tions, we suppress this effect. 

Another subtlety occurs for the square lattice system sizes 
L = 'in — 1, with n a positive integer number, for which 
the pattern of our pinning fields allows two different ordered 
states, corresponding to the two different orientation of the di- 
agonal stripes. 23 For these cases, at a sufficiently large bond 
dimension, a superposition of both types of order will occur 
and lead to a significant decrease of the local moment (ex- 
cept on sites where the two types of order coincide) and a 
significant increase of the entropy. In these cases, we pin two 
additional sites to uniquely select the order. 



4. Extrapolation 

The reliable extrapolation of results obtained for a lim- 
ited bond dimension to the limit of infinite bond dimension, 
where DMRG becomes exact, remains a challenge. While 
some reliable results are known for one-dimensional critical 
systems, ' 1 one has to resort to heuristic techniques in more 
general situations, such as two-dimensional systems. Such 
techniques include the extrapolation in the truncated weight 
or in the variance. Since we use a single-site optimization 
method, the truncated weight cannot be obtained reliably, and 
the calculation of the variance is only possible for smaller sys- 
tem sizes. We therefore resort to an extrapolation of the mag- 
netization in the bond dimension using the values obtained for 
the three largest values of D; for bond energies, we use only 
the result obtained for the largest value of D. While most 
simulations were performed using up to D = 4800 states, 
we have confirmed the accuracy of our results with up to 
D = 6400 states for some selected systems. 

Similarly, an accurate finite-size extrapolation is difficult 
given the few system sizes we can access. Also, the de- 
pendence of the order parameter and the energy on the sys- 
tem size, boundary conditions and aspect ratio is not known. 
In fact, previous studies have even observed surprising cases 
such as non-monotonic behavior for very small systems. 28 



C. Infinite projected entangled-pair states (iPEPS) 

1. Setup 

An iPEPS is a tensor network made of a set of rank-5 ten- 
sors periodically repeated on a two-dimensional lattice to effi- 
ciently represent ground state wave functions in the thermody- 
namic limit. ~ Each tensor has four auxiliary bonds which 
connect to the four nearest-neighbor tensors, and a fifth index 
carrying the local Hilbert space of a lattice site. The accuracy 
of the ansatz can be controlled by the dimension of the auxil- 
iary bonds, called the bond dimension D. As the optimization 
scheme for the tensors we perform an imaginary time evolu- 
tion with the so-called simple update (see Refs. 38 and 39) 
adopted from the time-evolving block decimation method in 



one dimension. 40,41 For the square lattice we verified the re- 
sults up to D = 8 also with the full update (see e.g. ), which 
is optimal but has a higher computational cost. The triangular 
lattice simulations are done with the same ansatz as for the 
square lattice, but now with an additional next-nearest neigh- 
bor interaction along one of the diagonal directions. The up- 
date scheme for this case is explained in Ref. 42. 

We performed simulations with different rectangular unit 
cells of size L x x L y in iPEPS . To represent the state with 
three-sublattice order efficiently, a 3 x 3 cell is used, with 
3 different tensors Ta,Tb, and Tq for the three sublattices 
respectively. We verified that the same state is obtained by 
using a similar cell with 9 different tensors. The 2x2 unit 
cell is used to enforce a state with two-sublattice order. 

To contract the tensor network efficiently, e.g. for the com- 
putation of observables, the corner transfer matrix scheme 
adapted to large unit cells 42 is used. The accuracy of the 
approximate contraction can be controlled by the so-called 
boundary dimension \- For large values of Z) a x up to 250 is 
used, where quantities of interest are extrapolated in x, with 
an extrapolation error being small compared to symbol sizes. 
For a better efficiency we use tensors with 7L q symmetry, a 
discrete abelian subgroup of SU(3). 45-47 



2. Calculation of the order parameter 

Since iPEPS is an ansatz for the wave function in the ther- 
modynamic limit, the SU(3) symmetry may be spontaneously 
broken, leading to a finite local moment m defined in Eq. 24. 
In order to pin the direction of the moment in SU(3) color 
space an initial field is applied, which is taken to zero at a 
later stage of the imaginary time evolution. We verified that 
we obtain the same results without initial field, and by com- 
puting the moment taking all generators of SU(3) into account 
(see Eq. (A4) in appendix A). 



3. Extrapolation 

For highly entangled systems quantities of interest such as 
the energy or the local moment are typically not converged 
as a function of the bond dimension D at the maximal value 
of D used, and thus an extrapolation to the infinite D limit 
is desirable. However, in general the dependence of observ- 
ables on D is (still) unknown, which limits the accuracy of 
such extrapolations. Since the approach is variational the en- 
ergy decreases with increasing D, and therefore the energy at 
the largest value of D provides an upper bound of the exact 
energy. Empirically, the exact value lies between the linear 
extrapolated value and the value at the largest D, and thus we 
take the middle of these two values as an estimate and the 
difference between the two values as an error bar. The same 
holds for the local moment, which is typically suppressed with 
increasing D, since more quantum fluctuations are taken into 
account with increasing D which renormalize the ordered mo- 
ment. Typically, the energy converges faster than the local 
moment. 
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FIG. 2. Comparison of the energy per site (upper panel) and the local 
moment (lower panel) of the SU(3) model (2) on the triangular lat- 
tice, obtained from LFWT, ED, DMRG and iPEPS. In each plot, the 
DMRG results are shown as a function of the inverse system length 
1/L (lower x-axis), whereas the iPEPS results are shown as a func- 
tion of inverse bond dimension 1/D (upper x-axis). Dotted lines are 
only guides to the eye. 



III. RESULTS FOR THE TRIANGULAR LATTICE 

We first present the results for the energy and the ordered 
moment for the model on the triangular lattice, obtained with 
LFWT, ED (energies only), DMRG and iPEPS. As mentioned 
in the introduction, on the triangular lattice a three sub-lattice 
order is already obtained from a simple product state ansatz. 
Inclusion of quantum fluctuations via LFWT does not destroy 
the order, 21 but renormalizes the local moment. In the fol- 
lowing we show that this holds even when including further 
quantum fluctuations. 




FIG. 3. Bond energies and local color densities in the triangular 
(7 + 2) x 7 lattice obtained from DMRG. The thickness of the bonds 
is proportional to the magnitude of the bond energy. An external po- 
tential is applied on the first and the last column to pin the sites to a 
specific color. 



E s « —0.69(1) J. 

Since we use open boundary conditions in the DMRG sim- 
ulations, the energy per site is not uniform in the system. Fig- 
ure 3 shows that the energy per bond close to the boundary 
is lower than far away from the boundaries. To obtain an es- 
timate of the energy per site in the "bulk" we average over 
the six bond energies around the central site for odd systems 
sizes. For even system sizes, we average over the four sites at 
the center of the system. This estimate is plotted in Fig. 2a) 
for different system sizes. The energy first increases with sys- 
tem size, and decreases slightly from L = 7 to L = 8. For the 
largest system L = 8 the estimated energy per site in the bulk 
isE s = -0.6775 J. 

Comparing the iPEPS energies obtained with different unit 
cell sizes, we find that the 3x3 unit cell yields a considerably 
lower variational energy than the 2x2 unit cell, which indi- 
cates that the symmetry breaking in the ground state is com- 
patible with the 3-sublattice order. The energy per site has 
not converged yet as a function of bond dimension D. Since 
the energy typically converges faster than linearly in 1/D we 
(empirically) expect the energy to lie in between the value for 
the largest D, E® =1 ° = —0.672 J, and the energy obtained 
from linear extrapolation of the last three data points in 1/D, 
El x = —0.708 J. Taking the mean of these two values yields 
an estimate of E s = —0.69(2), which is compatible with the 
DMRG result for the largest system. 



B. Local moment 



A. Energy per site 

Figure 2a) shows a comparison of the energy per site ob- 
tained with the four methods, where linear flavor-wave the- 
ory predicts a value of E s = — 0.6295J. Extrapolating 
ED energies for symmetric clusters consisting of N = 9, 
12, and 21 sites using a standard 1/N 3 / 2 form, we obtain 



In Fig. 2b) we present the results for the local moment m 
obtained with the various approaches, where to = 1 for the 
fully polarized case. 

As mentioned in Sec. II A, linear flavor-wave theory pre- 
dicts a value of to = 0.484. 

The DMRG results correspond to the local moment at the 
central site of the system for odd system sizes. For even sys- 
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tern sizes, the magnitude of the ordered moment is averaged 
over the four sites that make up the central plaquette of the 
system. Variations depending on the distance to the bound- 
aries in x- and y- direction can be observed, as shown in 
Fig. 4a). The value is decreasing with increasing distance 
from the pinning sites (in x-direction), whereas the value 
is seen to increase away from the boundary in y-direction. 
As a function of system size the local moment is increas- 
ing. As mentioned before, an accurate extrapolation to the 
thermodynamic limit is challenging, but a value in the range 
to s» 0.43 — 0.6 seems compatible with the DMRG data. 

The local moment obtained with iPEPS decreases with in- 
creasing D, an effect which can also be observed e.g. in the 
SU(2) Heisenberg model. With increasing D more quantum 
fluctuations are taken into account which reduce the magnetic 
moment from its value in the classical (product-state) limit, 
corresponding to D = 1. For the largest bond dimension used, 
D = 12, we find a value m — 0.58, whereas a linear extrap- 
olation in 1/D suggests a value of m = 0.52. As discussed 
in Sec. II C 3 the exact scaling behavior of m with 1/D is not 
known, but empirically, we expect m to lie in between these 
two values. 



C. Discussion 

Even though we cannot determine 777 up to a high preci- 
sion all three methods clearly suggest that the ground state 
has three-sublattice order with a large local magnetic moment 
in the range m = 0.43 — 0.6. As in the case of the SU(2) 
Heisenberg model on the square lattice, linear flavor wave the- 
ory (spin wave theory) already gives a good estimate of the 
local moment. 
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IV. RESULTS FOR THE SQUARE LATTICE 

We next consider the SU(3) Heisenberg model on the 
square lattice. As explained in the introduction a site factor- 
ized ansatz leads to an infinite number of degenerate ground 
states and quantum fluctuations (with LFWT) selects the 
three-sublattice state. Thus, quantum fluctuations seem to 
play a more dominant role on the square lattice, and it is con- 
ceivable that another ground state is selected when further 
quantum fluctuations beyond LFWT are taken into account. 
However, we show in the following that this is not the case 
here, i.e. that the three-sublattice order is stable and that addi- 
tional quantum fluctuations only further renormalize the local 
moment. 



A. Energy per site 

In Fig. 5a), the value of the energy per site from linear 
flavor-wave theory, E s = — 0.725J, has previously been 
calculated in Ref. 23, and is low compared to the numeri- 
cal results. We note the LFWT energy is not variational, so 
that it can be lower than the exact ground state value. We 



FIG. 4. DMRG results with D = 4800 states: Local moments as 
defined in Eqn. (24) for the triangular (top panel) and square (bottom 
panel) lattice for a system size (7 + 2) x 7. The plateau is very flat 
in the case of the square lattice, while corrections from the boundary 
are more pronounced for the triangular lattice. 



have also included an ED estimate of the energy per site 
E s = —0.63185 J, which is based on an extrapolation using 
square samples with N = 9 and 18 sites. 

The DMRG energy per site is E s = —0.625 for the largest 
system, and seems to further increase as a function of system 
size. As in the triangular lattice case we estimate the bulk en- 
ergy by taking the mean value over the bonds adjacent to a 
central site for odd system sizes and four sites for even system 
sizes. This energy seems higher than the LFWT and iPEPS 
result, which could indicate that boundary effects are large so 
that we do not get a good estimate for the "bulk" energy, or it 
could be that for larger systems the energy as a function of sys- 
tem size decreases again. We further note that an anisotropy 
in the bond energies can be observed, with stronger bonds in 
y-direction than in x-direction, shown in Fig. 6. 

Comparing the energies from different unit cell sizes in 
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FIG. 5. Same plot as in Fig. 2 but for the SU(3) model (2) on the 
square lattice. 



iPEPS, we observe a similar behavior as on the triangular lat- 
tice, namely that the 3x3 unit cell provides a better variational 
energy than the 2x2 unit cell for all values of D. The esti- 
mated energy per site in the limit D — >• oo is E s = —0.66(1). 




FIG. 6. Bond energies and local color densities in the square 
(7 + 2) x 7 lattice obtained from DMRG. The thickness of the bonds 
is proportional to the magnitude of the bond energy. An external po- 
tential is applied on the first and the last column to pin the sites to a 
specific color. 



in the limit D — > oo. The value for the largest bond dimen- 
sion, D = 16, is m = 0.3422 which is close to the DMRG re- 
sult for the largest system. However, in the limit D —> oo the 
data suggests a lower value of roughly m — 0.25(5), which is 
lower than the prediction from DMRG. 



C. Discussion 

Both the DMRG and iPEPS results are compatible with the 
proposed 3-sublattice Neel ordered ground state. From the 
present data we can only give a rough estimate of the ordered 
moment in the thermodynamic limit of m = 0.2 — 0.4, which 
is clearly finite, but smaller than on the triangular lattice. 



V. SUMMARY 



B. Local moment 

Figure 5b) summarizes our results for the local moment, 
obtained from DMRG and iPEPS. As explained in Ref. 23, 
the ordered moment cannot be computed within LFWT. 

The finite size effects observed in DMRG are qualitatively 
different from the triangular lattice case. The local moment as 
a function of distance of x in Fig. 4 reaches a plateau already 
after 3 sites away from the border, which could suggest that 
finite size effects on the ordered moment are smaller than on 
the triangular lattice. In Fig. 5b) the local moment in the mid- 
dle of the system first decreases and then increases with sys- 
tem size with, however, a smaller slope than in the triangular 
lattice case. Thus, the present data is compatible with a non- 
vanishing local moment in the thermodynamic limit, which is 
smaller than on the triangular lattice. The ordered moment of 
the largest system is m = 0.368. 

The local moment obtained with iPEPS decreases with in- 
creasing bond dimension but is not seen to extrapolate to zero 



Our study confirms that the ground state of the SU(3) 
Heisenberg model exhibits a three-sublattice order on both the 
triangular and the square lattice, in accordance with previous 
predictions by LFWT and exact diagonalization. 20 ' 21,23 

The situation on the triangular lattice resembles the one of 
the SU(2) Heisenberg model on the square lattice. In both 
cases the ground state can already be understood at the clas- 
sical level, and quantum fluctuations simply renormalize the 
ordered moment. These fluctuations are well captured already 
within linear flavor wave theory (i.e. spin wave theory in the 
SU(2) case). With iPEPS the ordered moment decreases with 
increasing bond dimension D, which can intuitively be un- 
derstood because the bond dimension controls the amount of 
quantum fluctuations taken into account. All three methods 
used in this study yield a finite ordered moment in the range 
m = 0.43 — 0.6. The uncertainty in this value stems from the 
error in the extrapolation to the thermodynamic limit in the 
case of DMRG, and from the extrapolation to the infinite D 
limit in the case of iPEPS. 
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In the case of the square lattice, the order cannot be pre- 
dicted at the classical level. Quantum fluctuations hence play 
a very different role than in the case of the triangular lat- 
tice: instead of renormalizing the mean-field result, they stabi- 
lize the three-sublattice order against other competing states. 
Quantum effects are therefore more important both qualita- 
tively and quantitatively, and an estimate of the ordered mo- 
ment in the thermodynamic limit has previously been lacking. 
Both DMRG and iPEPS predict a finite value in the range 
to = 0.2 — 0.4, i.e. the ordered moment is more strongly 
suppressed than on the triangular lattice, but clearly finite. 
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Appendix A: Schwinger boson study 

The Schwinger boson mean-field theory (SBMFT), intro- 
duced by Arovas and Auerbach and extended by Read 
and Sachdev, has been widely used for SU(2) models and 
more recently for models with SU(4) symmetry and SU(A^) 
models. 54,55 This approach is justified in the context of a large 
Af expansion, where Af is the number of boson flavors. Here, 
we stress that N and Af are two different numbers and the 
mean-field approach is equivalent to taking the limit Af to in- 
finity with N fixed. After a brief summary of the SBMFT, we 
address the following question: Is this theory really adapted 
to the study of models with SU(N) symmetry when N > 21 

We use the Schwinger bosons defined in Sec. II A and im- 
pose the constraint on the boson number at each lattice site i: 



(Al) 



For the model of Eq.(2), k = 1, and the Hamiltonian is 
the same as Eq. 7. We define the link operator Aij a p = 
(b a (i)bp(j) — b a {j)bp{i))/^f2 such that the permutation op- 
erator writes 



Pi 



1 



2 AU/ ; 

a>/3 



Aijot/3 • 



(A2) 



The Hamiltonian is thus of degree four in bosonic operators 
and is not directly solvable. A mean-field (MF) approximation 
lowers the degree to two: the Hamiltonian becomes quadratic 
and solvable by a Bogoliubov transformation. One possible 
MF parameter is A.y Q( g = (Aij a p). It is the most often used 



in SU(2) SBMFT because it is invariant by SU(2) transforma- 
tions : Aijap is the destruction operator of a SU(2) singlet of 
colors a and f3. The MF Hamiltonian writes 

#MF = 1 - 2 I] ( A *^l Q/ 3 + h.C. ~ l A ^l 

(ij) \ a>fi 

+ £)A(i)(K-n(i)), (A3) 

i 

where a Lagrange multiplier X(i) is used to impose the con- 
straint of Eq. Al on average at each site. 

The success encountered by this theory for SU(2) mod- 
els comes from the fact that two phases are possible for the 
ground state of the Hamiltonian of Eq. A3. For n lower 
than a critical value k c , the ground state is invariant by 
SU(2) global spin transformations and the elementary exci- 
tations are gapped spinons. If the Hamiltonian symmetries 
are respected, 56 ' 37 this phase is a topological spin liquid. For 
k > k c , one or several spinons become gapless, and to sat- 
isfy the constraint, a Bose condensate is required, breaking 
the SU(2) symmetry of the ground state: we have a long range 
ordered state. Unlike spin-wave expansions, this theory does 
not assume any order. The system is free to order or not, and 
the pattern is not imposed. 

This property comes from the fact that the MF Hamiltonian 
(Eq. A3) is expressed in terms of SU(2) invariant link opera- 
tors. For N > 2, singlets occupy TV sites and can no longer 
be destroyed by quadratic operators. This implies that the MF 
Hamiltonian is no longer SU(AT) invariant. This can be veri- 
fied using the order parameter to obtained using the Casimir 
operator on a site 



\ 



N-l 



(A4) 



o,j8 



For a product (non entangled) state, to = 1 and for a site in 
a SU(iV) singlet, to = . For SU(2), we recover to = |(S)|. 
For N > 2, the only way to have to = would be to set all 
the MF parameters Aij a p to 0. Thus even without any gapless 
spinon, the ground state is already long range ordered and we 
cannot have a spin liquid. The condensation is then only the 
breaking of a remaining freedom of the bosons. 

We now come back to the SU(3) model on the square and 
triangular lattice. The n — > oo limit is the classical limit, 
namely the 3 state Potts model. We first calculate the MF pa- 
rameters in this limit for the orders considered in this article, 
with 2 (or 3) sublattices denoted A, B (and C). These or- 
ders still have the degeneracies associated to the global SU(3) 
symmetry. Depending on the choice of 3 orthogonal vectors 
Ux, X = A, B,C, which specify the orientation on the dif- 
ferent sublattices, the mean-field parameters Aij a p take dif- 
ferent phases and modulus, unlike SU(2) where the only de- 
gree of freedom was the gauge choice (the values of Aij a p did 
not depend on how SU(2) was broken). We once again note 
that a SU(3) spin liquid cannot be described in this formalism. 



10 



lattice 


state 


-Bmf 


m 


n c 


triangular 


three sublattices 


-0.58 


0.93 


0.86 


square 


two sublattices 


-0.68 


0.50 


0.61 


square 


three sublattices 


-0.62 


0.89 


0.73 



TABLE I. Results obtained by SBMFT. Emf is the energy per site 
obtained by Eq. A3 with mean-field parameters verifying the self- 
consistency conditions, m is defined in Eq. A4 and n c is the number 
of bosons in the gapless modes(s). These quantities are extrapolated 
to the thermodynamical limit. 



Let us choose 



In the k — V oo limit, we can replace the b a (i) operators 
by their mean values (b a (i)}. If i is on the X sublattice, 
{b a (i)) — (ux)a- Thus we obtain the MF parameters 
for all the considered orders: Ay Q( 9 = ((b a (i))(b/3(j)) — 

(b a (j))(b,m/V2. 

For finite n, these parameters are not self-consistent and 
have to be adjusted. The chemical potential A is assumed to 
be site independant. We restrict our search for mean-field so- 
lutions to states obtained from the classical states by changing 
the modulus of the non zero A's. The energy, the order pa- 
rameter and the fraction of condensed bosons n c are given in 
Tab. A for the states discussed in this article. In all cases, m is 
far larger than the order parameter obtained in LFWT, DMRG 
and iPEPS, even without any boson condensation. The ener- 
gies obtained are not variational since the boson Hilbert space 
is larger than the physical one (the ground state is a superpo- 
sition of states with different boson numbers on each site). 

The three sublattice state on the square lattice has a larger 
energy than the two sublattice state, in contrast to the results 
obtained with LFWT, ED, DMRG and iPEPS. This result can 
be understood qualitatively. The SU(3) two-sublattice and the 
SU(2) two-sublattice state (historically treated by Arovas and 
Auerbach ) share several properties: same energy and same 
condensed fraction. The value of to = 0.5 of the order pa- 



rameter for SU(3) is exactly the value for a site in a SU(2) 
singlet and corresponds to m = for SU(2). By fixing the 
MF parameters, a direction among the three available is for- 
bidden to the bosons, that are confined in a SU(2) manifold. 
This remaining symmetry can only be broken by a conden- 
sate. The SU(3) ground state is thus exactly the same as the 
SU(2) one. The only difference stays in the existence of addi- 
tional excitations in the third direction for SU(3). They have a 
maximal energy cost, the value of the chemical potential, and 
form a flat band over all other excitations. Thus, the energy 
cannot be lowered by fluctuations in the full SU(3) space. We 
see in Tab. A that the magnetization is even larger in the other 
states (the three sublattice states on the square and triangular 
lattice), with to = 0.89 and 0.93 as compared to to = 0.5. 
The fluctuations are even more constrained, which provides a 
plausible explanation why the energy of the three-sublattice 
state on the square lattice is larger than the energy of the two 
sublattice state. 

To conclude this appendix, let us put these results in a 
broader perspective. For SU(2), the SBMFT is a convenient 
way to go beyond linear spin-wave theory (LSWT). This for 
instance allows to lift the classical degeneracies that may sur- 
vive LFWT. The kagome antiferromagnet offers a good ex- 
ample. Within LSWT, coplanar classical ground states remain 
degenerated, whereas SBMFT lifts the degeneracy in favor of 
the \/3x For that model, a classical state with higher lin- 
ear spin-wave energy has recently been shown to have an even 
lower SBMFT energy. In the context of the SU(3) model on 
the square lattice, the motivation to use SBMFT was to re- 
move the line of soft modes obtained in LFWT for the three- 
sublattice order, and at the same time the divergence of the 
correction to the magnetization. This is indeed achieved by 
the SBMFT, but unfortunately the resulting picture is not sat- 
isfactory: For SU(iV) models with A^ > 2, a spontaneous 
breaking of the SU(A0 symmetry is not possible since the 
choice of MF parameters already breaks it. No spin liquid 
ground state exists for the MF Hamiltonian of Eq. A3 and 
the quantum fluctuations taken into account in a long-range 
ordered ground-state are limited to some subspace of SU(3). 
As a consequence, we cannot use it to compare the MF states 
derived from several SU(3) classical orders. Other bosonic 
representations of SU( AO spins could lead to better results. 
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